library(foreign)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(stargazer)
library(pscl)
library(ggeffects)

#laptop
setwd("/Users/joelsievert/Dropbox/Research/Legislative Responsiveness/National Road/data/PRQ Replication")

#desktop
#setwd("/Users/jsievert/Dropbox/Research/Legislative Responsiveness/National Road/data/roll call votes")

####################################
##votes to extend road##
###################################

h21 <- read.csv("h21_extend.csv")
h23 <- read.csv("h23_extend.csv")

#national road

m1 <- glm(rc_natl ~ jackson  + dist_road + dist_natl + slave_p +  distrib_cmte, family = binomial, data = h21)
m2 <- glm(rc_natl ~  dist_road + dist_natl + slave_p + distrib_cmte , family = binomial, data = subset(h21, jackson == 1))

o1 <- coeftest(m1, vcov = vcovHC(m1, "HC0"))
o2 <- coeftest(m2, vcov = vcovHC(m2, "HC0"))

p1 <- ggpredict(m1, terms = "jackson", condition = c(dist_road = 2.1, dist_natl = 3.3, slave_p = 0, distrib_cmte = 0))
p2 <- ggpredict(m2, terms = "distrib_cmte", condition = c(dist_road = 2.1, dist_natl = 3.3, slave_p = 0) )


#maysville
m3 <- glm(rc_mays ~ jackson + dist_road + dist_mays + slave_p +  distrib_cmte, family = binomial, data = h21)
m4 <- glm(rc_mays ~ dist_road +  dist_mays + slave_p +  distrib_cmte , family = binomial, data = subset(h21, jackson == 1))

o3 <- coeftest(m3, vcov = vcovHC(m3, "HC0")) 
o4 <- coeftest(m4, vcov = vcovHC(m4, "HC0")) 

p1 <- ggpredict(m3, terms = "jackson", condition = c(dist_road = 2.1, dist_natl = 3.3, slave_p = 0, distrib_cmte = 0) )
p2 <- ggpredict(m4, terms = "distrib_cmte", condition = c(dist_road = 2.1, dist_natl = 3.3, slave_p = 0) )


##23rd
m5 <- glm(rc ~ jackson   + dist_road + dist_natl + slave_p + distrib_cmte , family = binomial, data = h23)
m6 <- glm(rc ~ dist_road + dist_natl + slave_p + distrib_cmte , family = binomial, data = subset(h23, jackson == 1))

o5 <- coeftest(m5, vcov = vcovHC(m5, "HC0")) 
o6 <- coeftest(m6, vcov = vcovHC(m6, "HC0")) 

p1 <- ggpredict(m5, terms = "jackson", condition = c(dist_road = 2.38, dist_natl = 9.6, slave_p = 0, distrib_cmte = 0) )
p2 <- ggpredict(m6, terms = "distrib_cmte", condition = c(dist_road = 2.38, dist_natl = 9.6, slave_p = 0) )


stargazer(o1, o2, o3, o4, o5, o6, digits = 2)


#calculate PRE
#See code from Dave Armstrong
pre(m1)
pre(m2)
pre(m3)
pre(m4)
pre(m5)
pre(m6)






